##################################################################
# File:      model_fitting_online_appendix.R
# Purpose:  Replication of the online appendix figures and tables for 
#                 Journal of Conflict Resolution
# Author:    Nadiya Kostyuk
# R version 4.2.2 (2022-10-31 ucrt)
# Last Edit: 01/05/2023
##################################################################


rm(list=ls())

## Install & load packages (all at once)
list.of.packages <- c('knitr', 'countrycode', 'dplyr','foreign','glmmTMB','lme4',
                      'stargazer')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]; if(length(new.packages)){install.packages(new.packages,dependencies=TRUE)}
lapply(list.of.packages, require, character.only = TRUE)

# setting directory: 
setwd("")


# loading the data:
load('main_data.RData')


################################################
# Selecting the model that provides the best fit
# (Section 3 of the Online Appendix)
# Table 2
# #############################################



# creating an empty list to save results: 
mod_list = list()

# 1. no random effects 
mod_st1 <- glm(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA),
               data=data,
               family=binomial)
summary(mod_st1)
mod_list[[1]] = mod_st1


# 2. accounting for dependence within dyads (does not account for type):
mod_st2 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA)
                       +(1|DYAD),
                       data=data,
                       family=binomial,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod_st2)
mod_list[[2]] = mod_st2


# 3. accounting for dependence within dyads for each type of conflict: THIS IS THE MAIN MODEL used in the manuscript
mod_st3 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA)
                       +(1+TYPE2||DYAD),
                       data=data,
                       family=binomial,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod_st3)
mod_list[[3]] = mod_st3


# 4. accounting for dependence within dyads for each type of conflict, allows for positive correlation within dyad-year
mod_st4 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA)
                       +(1+TYPE2||DYAD)+(1|DYAD_YEAR),
                       data=data,
                       family=binomial,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod_st4)
mod_list[[4]] = mod_st4


# 5. accounting for dependence within dyads for each type of conflict, allows for negative correlation within dyad-year
mod_st5 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA)
                       +(1+TYPE2||DYAD)+(0+TYPE2|DYAD_YEAR),
                       data=data, 
                       family=binomial,
                       control=glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e6))
)
summary(mod_st5)
mod_list[[5]] = mod_st5

# 6. accounting for dependence within dyads for each type of conflict, allows for positive and negative correlation within dyad-year
mod_st6 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA)
                       +(1+TYPE2||DYAD)+(1+TYPE2||DYAD_YEAR),
                       data=data,
                       family=binomial,
                       control=glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e6))
)
summary(mod_st6)
mod_list[[6]] = mod_st6

# saving model results
#save(mod_list, file = 'mod_list_tab2a.RData')





#########################
# Main Results (Table 3)
#########################

# 1. adding covariates
mod_st2a <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA)
                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                        +CINC_A_sc+INV_CAP_DISTANCE_sc
                        +NUCLEAR
                        +(1+TYPE2||DYAD),
                        data=data,
                        family=binomial,
                        control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod_st2a)
mod_list2 = list()
mod_list2[[1]] = mod_st2a


# 2. adding interactions with covariates (i.e., internet users per capita for the attacker and the target)
mod_st2b <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc)
                          +CINC_A_sc+INV_CAP_DISTANCE_sc
                          +NUCLEAR
                          +(1+TYPE2||DYAD),
                          data=data,
                          family=binomial,
                          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod_st2b)
mod_list2[[2]] = mod_st2b


# 3. adding interactions with covariates (i.e., internet users per capita for the attacker and the target, cinc score for the attacker)
mod_st2c <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                        +CINC_A_sc)+INV_CAP_DISTANCE_sc
                          +NUCLEAR
                          +(1+TYPE2||DYAD),
                          data=data,
                          family=binomial,
                          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod_st2c)
mod_list2[[3]] = mod_st2c


# 4. adding interactions with covariates (i.e., internet users per capita for the attacker and the target, cinc score for the attacker, distance between capitals)
mod_st2d <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                        +CINC_A_sc+INV_CAP_DISTANCE_sc)
                          +NUCLEAR
                          +(1+TYPE2||DYAD),
                          data=data,
                          family=binomial,
                          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod_st2d)
mod_list2[[4]] = mod_st2d

# 5. adding interactions with all covariates 
mod_st2e <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                        +CINC_A_sc+INV_CAP_DISTANCE_sc
                                        +NUCLEAR)
                          +(1+TYPE2||DYAD),
                          data=data,
                          family=binomial,
                          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod_st2e)
mod_list2[[5]] = mod_st2e


############################################
# Robustness Checks
###########################################

###########################################
# Section 4.1 of the Appendix
# See: reporting_biases_online_appendix.R
############################################

########################
# 4.2 (Model)
data=loadRData('main_data.RData')
mod1 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                 +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                 +CINC_A_sc)+INV_CAP_DISTANCE_sc
                   +NUCLEAR
                   +(1+TYPE2||DYAD),
                   data=data,
                   family=binomial,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod1)

mod2 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                  +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                  +CINC_A_sc)+INV_CAP_DISTANCE_sc
                    +NUCLEAR
                    +(1+TYPE2||DYAD)+(1|DYAD_YEAR),
                    data=data,
                    family=binomial,
                    control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod2)

mod3 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                  +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                  +CINC_A_sc)+INV_CAP_DISTANCE_sc
                    +NUCLEAR
                    +(1+TYPE2||DYAD)+(0+TYPE2|DYAD_YEAR),
                    data=data,
                    family=binomial,
                    control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod3)

mod4 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                  +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                  +CINC_A_sc)+INV_CAP_DISTANCE_sc
                    +NUCLEAR
                    +(1+TYPE2||DYAD)+(1+TYPE2||DYAD_YEAR),
                    data=data,
                    family=binomial,
                    control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod4)



########################
# 4.3 Month
load('data_month.RData')
mod <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                        +CINC_A_sc)+INV_CAP_DISTANCE_sc
                          +NUCLEAR
                          +(1+TYPE2||DYAD),
                          data=data,
                          family=binomial,
                          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod)

########################
# 4.4 Distance
data=loadRData('main_data.RData')
mod <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                      +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                      +CINC_A_sc)+CONTIGUITY_4
                        +NUCLEAR
                        +(1+TYPE2||DYAD),
                        data=data_st,
                        family=binomial,
                        control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod)

########################
# 4.5 major powers:

mod <- lme4::glmer(OUTCOME~TYPE*((CYBER_LAG_NA+CONVENTIONAL_LAG_NA)*major_a
                                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                        +CINC_A_sc)+INV_CAP_DISTANCE_sc
                          +NUCLEAR
                          +(1+TYPE2||DYAD),
                          data=data,
                          family=binomial,
                          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod)

######################
# 4.6. Lags
####################

# reading the data: 
load('main_data_2dlag.RData')
mod <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                 +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                 +CINC_A_sc)+INV_CAP_DISTANCE_sc
                   +NUCLEAR
                   +(1+TYPE2||DYAD),
                   data=data,
                   family=binomial,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod)

##############################
# 4.7. Types of Operations
##############################


load('data_disrupt_deg.RData')
mod <- lme4::glmer(OUTCOME~TYPE*(DISRUPT_DEGRAD_LAG_NA+CONVENTIONAL_LAG_NA
                                 +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                 +CINC_A_sc)+INV_CAP_DISTANCE_sc
                   +NUCLEAR
                   +(1+TYPE2||DYAD),
                   data=data_disrupt_deg,
                   family=binomial,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod)


##############################
# 4.8. Additional Controls
##############################
mod <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                 +LOG_GDP_PERCAP_A_sc+LOG_GDP_PERCAP_B_sc
                                 +CINC_A_sc)+INV_CAP_DISTANCE_sc
                   +NUCLEAR
                   +(1+TYPE2||DYAD),
                   data=data,
                   family=binomial,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod)


mod <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                 +LOG_GDP_PERCAP_A_sc+LOG_GDP_PERCAP_B_sc
                                 +CINC_A_sc)+INV_CAP_DISTANCE_sc
                   +NUCLEAR
                   +(1+TYPE2||DYAD),
                   data=data,
                   family=binomial,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod)



